The 28th International Cosmic Ray Conference 



1 



Numerical Likelihood Analysis of Cosmic Ray Anisotropies 



m 
O 
O 

(N 



Carlos Ho j vat, 1 Thomas P. McCauley, 2 Stephen Reucroft, 2 and John D. Swain 2 

(1) Fermi National Accelerator Laboratory, P.O. Box 500, Batavia, IL 60510 

(2) Department of Physics, Northeastern University, Boston, MA 02115 



Abstract 



> 

o 

(N 

in 
o 

CO 

O 

6 

u 

in 
c3 



A numerical likelihood approach to the determination of cosmic ray anisotropies 
is presented which offers many advantages over other approaches. It allows a wide 
range of statistically meaningful hypotheses to be compared even when full sky 
coverage is unavailable, can be readily extended in order to include measurement 
errors, and makes maximum unbiased use of all available information. 

1. Introduction 

The search for anisotropies in the cosmic ray distribution (and of other 
objects - see for example [4,5]) is of increasing interest with the advent of the 
Pierre Auger Observatory (for a review of the physics up to this time, see [1]). 
The issues involved are, however, subtle and complicated by limited statistics at 
the highest energies and nonuniform sky coverage. The first full-sky search for 
UHECR anisotropies using a standard approach is presented in this conference [2] . 
This paper presents a maximum likelihood approach which is well-adapted to 
further studies in anticipation of much larger data sets. 

2. Expansion in Spherical Harmonics 



It is convenient to define a set of real spherical harmonics {^,m} f° r 
0,1,2,... andm = -£, -(£ - 1), 0, ...,(£- 1), I by 



A;| m 'pJ m '(cos#) cos(m0) for m = —£, . . . , —1, 
fc°P n °(cosfl) form = 0, 

kf P7 1 (cos 6) sm(m<b) for m = !,...,£ 



which are orthonormal with respect to the usual measure sm(9)d9d(fi and inte- 
gration over 9 from to tt and from to 2ir, the P™ are associated Legendre 
polynomials of the first kind, and the normalization constants k l m are: 
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The natural measure of anisotropy for a spherical distribution is in terms of these 
spherical harmonics, as each i labels, in a coordinate-independent fashion, just 
how much of each irreducible 50(3) representation is present. 

In a perfect world with infinite statistics and complete sky coverage there 
are now many possible approaches to estimating how much of each of these com- 
ponents is present in a distribution, or, better, what is the likelihood that a given 
function f(9,4>) with Fourier-Legendre expansion f(9,(p) = Efco^m=-< a «,»i^,m 
representing the probability density of sources gives rise to the observed distribu- 
tion g(9, 0). The coefficients can be extracted from the usual integral 

a t , m = r f(9, <j>)il>i, m (e, 0) sin Qd6d4 (3) 
Jo Jo 

Of course this gives no measure of what sort of error should be associated with 
the determined values of each coefficient. Alternative approaches are to fit for the 
coefficients by minimizing some x 2 -like quantity like 

f 2n «r g(9,<t>) -T,Zo^m=-e a £,mAm( 9 ^) . njnjl 

/ / snxddddd) (4) 

Jo Jo a 1 

with a a suitable measure of error, or to compute and maximize a corresponding 
likelihood that the hypothesized distribution parametrized by the a^ m gives rise 
to the the observed distribution g. Of course this is statistically unreasonable for 
a finite number of sources (i.e. to provide an infinite number of coefficients!), In 
general a decision must be made to truncate the expansion at some value of £ 
to make the sum finite, but now the general phenomenon of aliasing risks that, 
for example, a fit allowing for I = 0, 1 might give misleading results for observed 
data which is drawn from a purely 1 = 1 distribution, say. 

A more serious problem is that should part of the sky be unobserved, there 
is now no way to calculate anything! This is not a trivial point. An attempt to find 
anisotropies based only on observations in the Northern hemisphere with zeroes 
inserted for the whole Southern hemisphere would be wildly in error if some simple 
extrapolation were made to the unobserved region of the sky - especially if there 
were something bright and as-yet undetected in the South! The real challenge is 
to say something statistically meaningful with the data that is actually available. 
Clearly, observing only the Northern hemisphere and seeing a good degree of 
isotropy should increase one's net belief in overall isotropy of the full sky, while 
leaving open the possibility of a staggeringly bright or empty sky in the South. 
One approach is to try to make functions which would be orthonormal over the 
observed region of the sky, but it's not clear that this has much physical meaning 
as it elevates lack of acceptance to a status comparable to the SO (3) invariance of 
space embodied in the ipe,m- The following is a proposal for what seems to make 
good statistical sense. 
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3. A Likelihood Proposal 

Based on the above observations, we make the following proposal: keep the 
spherical harmonics as always with the (necessarily) truncated Fourier- Legendre 
expansion and construct an unbinned likelihood [3] function in which one clearly 
specifies which values of i are included in the sum. An unbinned likelihood func- 
tion, as described below, automatically makes maximum use of all detected in- 
formation, allows for measurement errors to be included easily, and is easy to 
implement numerically. Most importantly, however, the likelihood is not to be 
normalized as it stands. Rather, we take the Bayesian approach to likelihood 
which says that likelihoods give us ways to update our prior experimental data or 
guesses in light of new information. This will mean that we are able to present re- 
sults on various hypotheses about data without bias, and with the easy inclusion 
of other data from the same, or other experiments. 

To be concrete, we specialize here to the case where g is a sum of delta 
functions representing sources % of unit intensity at 6{, (pi and later discuss how 
one can treat the case of these being at uncertain locations, or taking into account 
other properties such as intensity, energy, composition, etc.. These will appear as 
natural generalizations to the approach described. 

The (unnormalized!) likelihood L({£}|#;, </>;) that the measured 8i,<f)i arise 
from f(9,4>) = Y^TT^ m =-i a t-,rn^i,m where the sum over i is specified according 
to whatever hypothesis is being tested (i.e. just taking {£} = {0, 1} allows for 
uniform and dipole contributions and no others, while {£} = {2} would be pure 
quadrupole) is 



and over the range in which the parameters are allowed to vary. An important 
caveat is that for / to represent a sensible probability distribution it should never 
be negative, and this must be checked. Two approaches are possible: one is to 
restrict the domain over which coefficients range so that the function is strictly 
positive (this is not actually very difficult in practice since distributions are often 
nearly uniform with small fluctuations superimposed) or to take as a probability 
distribution function some positive function of / in place of / above. In astron- 
omy [4] it is not uncommon to see exp(/). The choice ultimately represents the 
unavoidable presence of some (often hidden) assumption about what a sensible 
prior (i.e. in the absence of data) is for the likelihood - a point to which we 
return later. 

By construction then L({£}\6i,<fri) is normalized so that its integral over 
all the parameters that can vary (the coefficients a^ m included in the truncated 
Fourier- Legende expansion, and the total likelihood Lj-ot, which is a function of 




(5) 



where the integral in the denominator is over all the parameters that can vary 
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those a^ m is 

w(m) = iwmi^) (6) 

% 

This is a relative likelihood, and while no absolute normalization is possible (nor 
should it be!) if part of the sky is unobserved, it is now very useful for two types of 
calculation. If one has a prior expectation for the distribution of the a^ m (which 
might be that they are all a priori equally likely) then Ltot can be multiplied by 
this and the product treated as a normalized likelihood distribution for the Of m 
themselves. This is, of course, potentially dangerous, but does allow one to see 
how the new data (the measured (#,, (pi)) should cause one to revise earlier beliefs. 
Such a likelihood can be maximized with respect to the a>i. m and if not enough 
data is present to test the hypothesis (for example, more coefficents to fit for than 
data points) the fit will respond by just not converging (that is to say, there will be 
flat directions in the likelihood as a function of the parameters meaning that one 
can't decide) - the beauty of this approach is that it is, by construction, correct. 
When results are obtained, one automatically gets the values of the parameters, 
and their entire likelihood distributions from which errors (which need not even 
be Gaussian) can be extracted. One can even do things like fix some parameters 
to those given by a favoured theory and then repeat the fitting and then obtain 
likelihoods for the correctness of that theory. 

More objectively, one can compute relative likelihoods in which the prior 
drops out, so that it is reasonable to ask (even in the absence of full sky coverage!) 
what the relative likelihood Lrel is of pure dipole distribution to one admitting 
uniform, dipole and quadrupole components: Lrel = Ltot ({0^2}) ^ one wan ^ s 
statements made to be about all energies over x, then one just uses the data points 
with energies over x. Similarly, data can be selected by composition and (even 
relative) anisotropies be searched for as functions of composition, energy, time, 
etc. Extensions to uncertainty in direction are trivial to include: simply divide a 
given event into a large number, M, say of subevents distributed appropriately and 
count each in the likelihood with a weight 1/M - this numerically convolves this 
uncertainty with the likelihood, and the size of any additional errors introduced 
by the procedure can be studied by varying M. 
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